Pair Trading with Machine Learning¶
Pair Trading - Introduction¶
Pair trading is a commonly used quantitative trading strategy used in hedge funds and investment banks. A major aspect of this strategy is the spread analysis/prediction amongst a pair of selected stocks.
The prinicple is as follows.
Let's say you have a pair of securities X and Y that have some underlying economic link. An example might be two companies that manufacture the same product, or two companies in one supply chain. If we can model this economic link with a mathematical model, we can make trades on it.
Process Flow¶
#Importing Dependencies
import numpy as np
import pandas as pd
import os
import glob
import pickle as pickle
import json
import missingno
import itertools
from itertools import cycle
from datetime import datetime
from builtins import object
from datetime import date
from dateutil.relativedelta import relativedelta
from collections import defaultdict
from concurrent import futures
import pandas_datareader.data as web
import warnings
warnings.filterwarnings('ignore')
#visualization dependencies
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
from plotly import graph_objects as go
from plotly import offline as py
import plotly.express as px
%matplotlib inline
#machine learning/statistical dependencies
import statsmodels.api as sm
from kneed import KneeLocator
from sklearn import metrics
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.metrics import mean_squared_error, silhouette_score
from sklearn.cluster import KMeans, AgglomerativeClustering, AffinityPropagation
from sklearn.manifold import TSNE
import statsmodels.api as sm
import statsmodels.tsa.stattools as ts
from statsmodels.tsa.arima_model import ARIMA
from pykalman import KalmanFilter
from keras.models import Sequential, Model
from keras.layers import Dense, LSTM, Dropout
import scipy.cluster.hierarchy as shc
Data Gathering¶
We scrape the web to gather the financial data of S&P 500 companies for last 5 years from current date to be utilized for the proposed pair trading strategy implementation.
Data Sources
The ticker, sector and industry metadata for S&P 500 companies was scraped from Wikipidea. (https://en.wikipedia.org/wiki/List_of_S%26P_500_companies)
The financial information for each company was gathered using python's
pandas_datareaderlibrary using the company tickers as identifiers.
We developed the class gather_data, which encapculates all the function utilized for the data gathering process.
Quick overview of the functions:
download_s_and_p_stats(): Used to fetch S&P 500 company metadata (ticker, sector, industry) from wikipidea.download_stock_financials(): Used to fetch financial data (open, low, high, close, volume, adj close) given a stock ticker.MultiThreader(): Used to speed up the process of fetching stock data using multi threading.combine_stock_data(): Used to merge all the stock CSV data files into on combined file.
The data gathering process consists of the following steps:
The
MultiThreader()is called through a class object. It utilizes thedownload_s_and_p_stats()to get the stock metadata. This metadata is then passed to thedownload_stock_financials()with 50 maximum parallel workers.The
download_stock_financials()extracts the ticker from each stock's metadata. This ticker is used with the start and end dates to fetch financial data of the stock for that time period.The financial data is merged with the metadata for each stock and a CSV file for each stock is created and stored at the given location.
We utilize the
combine_stock_data()to merge all the gathered stock CSV files into one dataframe and optionally store it in the given location
class gather_data(object):
def __init__(self):
super(gather_data, self).__init__()
#The start date for historical data
self.start_time = str(date.today()- relativedelta(years=5))
#The end date for historical data
self.end_time = str(date.today())
def download_s_and_p_stats(self):
'''
Fetches comapny ticker, sector and industry data for
S&P 500 companies
'''
data = pd.read_html('https://en.wikipedia.org/wiki/List_of_S%26P_500_companies')
df = data[0].reset_index()
#Extracting required features
df = df[['Symbol','GICS Sector','GICS Sub-Industry']]
#Formatting the group symbols
df['Symbol'] = [sym.replace('.','-') if sym.find('.') != -1 else sym for sym in list(df['Symbol']) ]
return df
def download_stock_financials(self,stock):
'''
Takes a stock ticker as input and fetches related financial data given the start and end time
Inputs
stock-- The stock ticker
The function generates and stores a csv file for stock financial data under the given path
'''
try:
stock_data = web.DataReader(stock[0],'yahoo', self.start_time, self.end_time)
stock_data[['Symbol','GICS Sector','GICS Sub-Industry']] = pd.DataFrame([stock], index=stock_data.index)
output_name = stock[0] + '_data.csv'
stock_data.to_csv(self.op_path +str(output_name))
except:
print('bad: %s' % (stock[0]))
def MultiThreader(self, op_path):
'''
We utilize the ThreadPoolExecutor from concurrent.futures module
to speed up the downloads utilizing parallelism
Input
op_path-- The path for output csv to be stored
'''
self.op_path = op_path
max_workers = 50
sp_stocks = self.download_s_and_p_stats()
#To handle case when smaller number of stock than thread are passed
workers = min(max_workers, len(sp_stocks))
with futures.ThreadPoolExecutor(workers) as executor:
res = executor.map(self.download_stock_financials, sp_stocks.values)
def combine_stock_data(self,curr_path, files_path, output_path = None, save_csv=False):
'''
Function to combine csv files for all stocks into one
Inputs:
curr_path-- Path for Current working directory
files_path-- Path where the files to be merged are stored
output_path-- Path to store the merged CSV file (default=None)(optional)
save_csv-- Boolean to specificy is CSV file should be stored (default=False)(optional)
'''
#changing the working directory to files path
os.chdir(files_path)
#fetching all file names with an '.csv' extension
file_extension = '.csv'
all_filenames = [i for i in glob.glob(f"*{file_extension}")]
#combining the files
combined_stocks_data = pd.concat([pd.read_csv(f ) for f in all_filenames])
if save_csv == True:
#changing working directory to the given output path
os.chdir(output_path)
#saving our combined csv data
combined_stocks_data.to_csv('combined_stocks_data.csv')
#chaning work directory back to original
os.chdir(curr_path)
return combined_stocks_data
if __name__ == '__main__':
#defining working directory paths
curr_path = os.getcwd()
files_path = curr_path + '/Data/Individual Stocks/'
output_path = curr_path + '/Data/'
#declaring gather_data class object
gd_obj = gather_data()
#downloading stock data
gd_obj.MultiThreader(files_path)
#combining the data files
combined_data = gd_obj.combine_stock_data(curr_path,files_path, output_path, True)
#displaying
combined_data.head()
#Shape of the combined data
np.shape(combined_data)
Data Preprocessing¶
The data preprocessing part contains of 2 steps:
Feature Engineering: We apply feature engineering on the data to extract multiple technical indicators which might help us to understand and realize the link between two stocks. For this we make use of
generate_technical_indicator()from thefeature_engineeringclass to derive and add aditional features to a given the stock financials data.Data Cleaning: The NA values are handles and the data is standardised to prepare it for future analysis
Technical Indicators¶
- RSI - Relative Strength Index, a momentum indicator used by traders to analyze bullish behaviou when a stock is oversold or bearish when a stock is overbought
$$\text{RS} = -\frac{AvgGain_{prev} * 13 + Current Gain}{AvgLoss_{prev} * 13 + Current Loss}$$
- KAMA - Kufman's Adaptive Moving Average, accounts for noise and volatile nature of the market. It closely follows price when noise is low and smooths out the noise when price fluctuates. It can be used to identify the turning points and price movements
- ADI - Accumulation/Distribution index acts as a leading indicator of price movements relative to the volume
$$\text{accdist} = accdist_{prev} + volume * CLV $$
- VPT - Volume Price trend is an indicator based on a running cumulative volume that adds or substracts a multiple of the percentage change in share price trend and current volume, depending upon the investment’s upward or downward movements.
- ATR - Average True Range, is an indicator that provides the intution of the degree of price volatility. $$\text{ATR} = \frac{1}{n}\sum_{i=1}^{n}TR_{i}$$ $$\text{TR} = max[(high-low), abs(high - close_{prev}),abs(close_{prev} - low)]$$
- Bollinger M-Band - Simple moving average of the price movements
- EMA - Exponential Moving Average of the price movements
Average Directional Movement Index - It is an indicator of strength of the trend for a given stock. . It is compose of two indicators, positive direction indicator (+DI) and negative direction indicator(-DI)
Moving Average Convergence Divergence - It is a trend-following momentum indicator that shows the relationship between two moving averages of prices in a given time window
Daily returns - The percentage gain from previous day's Adjacent closing price
Log returns - The percentage gain from logarithmic previous day's closing price
Reference: Wikipedia
class feature_engineering(object):
def __init__(self):
super(feature_engineering, self).__init__()
def generate_technical_indicator(self,data_df):
'''
Function to generate additional technical indicators for the stock
Input:
data_df-- Dataframe containing stock finacials data
Output:
Stock finacials data with added Dataframe of feature obtained from feature engineering
'''
# 1. Momentum Indicators
# Relative Strength Index
df = data_df
df['rsi'] = ta.momentum.rsi(df['Close'], window=14)
#Kaufman’s Adaptive Moving Average (KAMA)
df['kama'] = ta.momentum.kama(df['Close'],window=14)
# 2. Volume Indicators
# Accumulation/Distribution Index (ADI)
df['adi'] = ta.volume.acc_dist_index(df['High'], df['Low'], df['Close'], df['Volume'])
# Volume-price trend (VPT)
df['vpt'] = ta.volume.volume_price_trend(df['Close'], df['Volume'])
# 3. Volatility Indicators
# Average True Range (ATR)
df['atr'] = ta.volatility.average_true_range(df['High'], df['Low'],df['Close'], window=14)
# Bollinger Bands (BB) N-period simple moving average (MA)
df['bb_ma'] = ta.volatility.bollinger_mavg(df['Close'], window=20)
# 4. Trend Indicators
# Average Directional Movement Index (ADX)
df['adx'] = ta.trend.adx(df['High'], df['Low'], df['Close'], window=14)
# Exponential Moving Average
df['ema'] = ta.trend.ema_indicator(df['Close'], window=14)
# Moving Average Convergence Divergence (MACD)
df['macd'] = ta.trend.macd(df['Close'], window_fast=14, window_slow=30)
# 5. Other Indicators
# Daily Log Return (DLR)
df['dlr'] = ta.others.daily_log_return(df['Close'])
# Daily Returns
df['daily_returns'] = df['Adj Close'].pct_change()
# Moving Averages
averages = [50,200]
for avg in averages:
col_name = str(avg) +' Days Average'
df[col_name] = df['Adj Close'].rolling(window = avg, center = False).mean()
return df
#declaring feature_engineering class object
fe_obj = feature_engineering()
#Adding technical indicators
final_data = fe_obj.generate_technical_indicator(combined_data)
final_data.head()
#Null value check
final_data.isna().sum()
#Cleaning out the NA values
final_data = final_data.dropna()
np.shape(final_data)
#Null value check
final_data.isna().sum()
final_data.columns
#Applying Standard Scaling on the Contiunous features
train_data = final_data[['High', 'Low', 'Open', 'Close', 'Volume', 'Adj Close','rsi', 'kama', 'adi', 'vpt', 'atr',
'bb_ma', 'adx', 'ema', 'macd', 'dlr', 'daily_returns',
'50 Days Average', '200 Days Average']]
standard_df = StandardScaler().fit_transform(train_data)
standard_df.shape
#Null value check
train_data.isna().sum()
Pair Selection¶
Data Preparation for Pair Selection¶
Below is the data created after initial data scraping and feature engineering. It consists of 5 year daily data for all the S&P 500 tickers
final_data = pd.read_csv('/Users/prakharpatidar/Google Drive/DS/InfoViz DS 5500/Data/Work_data/final_data.csv', index_col = 0)
final_data = final_data.dropna().reset_index(drop=True)
final_data.head()
final_data.columns
df = final_data.pivot(index='Date', columns='Symbol', values='Adj Close')
df.head()
exp_data = df.copy()
pd.set_option('precision', 3)
exp_data.describe().T.head()
exp_data.isnull().values.any()
missingno.matrix(exp_data)
print('Data Shape before cleaning =', exp_data.shape)
missing_percentage = exp_data.isnull().mean().sort_values(ascending=False)
missing_percentage.head(10)
dropped_list = sorted(list(missing_percentage[missing_percentage > 0.2].index))
exp_data.drop(labels=dropped_list, axis=1, inplace=True)
print('Data Shape after cleaning =', exp_data.shape)
exp_data = exp_data.fillna(exp_data.median())
exp_data
missingno.matrix(exp_data)
returns = exp_data.pct_change().mean()*266
returns = pd.DataFrame(returns)
returns.columns = ['returns']
returns['volatility'] = exp_data.pct_change().std()*np.sqrt(266)
ret_data = returns
ret_data.head()
#Prepare the scaler
scale = StandardScaler().fit(ret_data)
#Fit the scaler
scaled_data = pd.DataFrame(scale.fit_transform(ret_data),columns = ret_data.columns, index = ret_data.index)
X = scaled_data
X.head()
Cluster Analysis¶
Cluster Analysis is a technique used to group sets of objects that share similar characteristics. It has been prominently used in the finance domain by investors to develop a cluster trading approach to help them build a diversified portfolio, helping protect the investor’s portfolio against systemic risks that could make the portfolio vulnerable to losses.
Here we explore three different clustering techniques for extracting different clusters of stocks exhibiting minimal correlation from one another. The three techniques are compared based on Silhouette Scores and the clusters from the technique with highest score is selected for further analysis. The pairs in each cluster are tested for cointegration via the ADF test giving us optimal pairs for our pair trading strategy.
K-Means Clustering¶
K-Means algorithm is an iterative algorithms that partitions the datasets into k distinct non-overlapping clusters where each datapoint belongs exactly to one cluster. It tries to cluster most similar (homogenous) points together while also maximizing distance between each individual cluster. The elbow method is used to iterate through the values of k and calculate the distortion for each value of k, and distortion and inertia for each value of k in the specified range, where distortion is the average of
the squared distances of each feature to the closest cluster center. We use Silhouette method to reaffirm the optimal number of clusters (K). The K-Means clustering model is trained using the found optimal K value to generate clusters for our pair selection.
K = range(1,15)
distortions = []
#Fit the method
for k in K:
kmeans = KMeans(n_clusters = k)
kmeans.fit(X)
distortions.append(kmeans.inertia_)
#Plot the results
fig = plt.figure(figsize= (15,5))
plt.plot(K, distortions, 'bx-')
plt.xlabel('Values of K')
plt.ylabel('Distortion')
plt.title('Elbow Method')
plt.grid(True)
plt.show()
kl = KneeLocator(K, distortions, curve="convex", direction="decreasing")
kl.elbow
#For the silhouette method k needs to start from 2
K = range(2,15)
silhouettes = []
#Fit the method
for k in K:
kmeans = KMeans(n_clusters=k, random_state=42, n_init=10, init='random')
kmeans.fit(X)
silhouettes.append(silhouette_score(X, kmeans.labels_))
#Plot the results
fig = plt.figure(figsize= (15,5))
plt.plot(K, silhouettes, 'bx-')
plt.xlabel('Values of K')
plt.ylabel('Silhouette score')
plt.title('Silhouette Method')
plt.grid(True)
plt.show()
kl = KneeLocator(K, silhouettes, curve="convex", direction="decreasing")
print('Suggested number of clusters: ', kl.elbow)
c = 8
#Fit the model
k_means = KMeans(n_clusters=c)
k_means.fit(X)
prediction = k_means.predict(X)
#Plot the results
centroids = k_means.cluster_centers_
fig = plt.figure(figsize = (18,10))
ax = fig.add_subplot(111)
scatter = ax.scatter(X.iloc[:,0],X.iloc[:,1], c=k_means.labels_, cmap="rainbow", label = X.index)
ax.set_title('k-Means Cluster Analysis Results')
ax.set_xlabel('Mean Return')
ax.set_ylabel('Volatility')
plt.colorbar(scatter)
plt.plot(centroids[:,0],centroids[:,1],'sg',markersize=10)
plt.show()
clustered_series = pd.Series(index=X.index, data=k_means.labels_.flatten())
clustered_series_all = pd.Series(index=X.index, data=k_means.labels_.flatten())
clustered_series = clustered_series[clustered_series != -1]
plt.figure(figsize=(12,8))
plt.barh(range(len(clustered_series.value_counts())),clustered_series.value_counts())
plt.title('Clusters')
plt.xlabel('Stocks per Cluster')
plt.ylabel('Cluster Number')
plt.show()
counts = clustered_series.value_counts()
cluster_vis_list = list(counts[(counts<20) & (counts>1)].index)[::-1]
clust1 = pd.DataFrame(clustered_series,columns=['label'])
grouped_df = clust1.groupby('label')
pivot_data = final_data.pivot(index='Symbol', columns='Date', values='Adj Close')
pivot_data.head()
#visualizing time series for clusters in K-means
for key, item in grouped_df:
print('Group: ',key)
grp_df = grouped_df.get_group(key)
Syms = grp_df.index
means = np.log(pivot_data.loc[Syms,].T.mean())
temp = np.log(pivot_data.loc[Syms,].T.sub(means))
# temp.plot()
# plt.show()
temp = temp.reset_index()
col_lst = temp.columns.values.tolist()
col_lst = col_lst.remove('Date')
df_melt = temp.melt(id_vars='Date', value_vars=col_lst)
fig = px.line(df_melt, x='Date', y='value',color='Symbol')
fig.show()